library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
-- Attaching packages ------------------------------------------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.5 v purrr 0.3.4
v tibble 3.1.2 v dplyr 1.0.7
v tidyr 1.1.3 v stringr 1.4.0
v readr 1.4.0 v forcats 0.5.1
-- Conflicts ---------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(spatstat)
Lade nötiges Paket: spatstat.data
Lade nötiges Paket: spatstat.geom
Registered S3 method overwritten by 'spatstat.geom':
method from
print.boxx cli
spatstat.geom 2.3-0
Lade nötiges Paket: spatstat.core
Lade nötiges Paket: nlme
Attache Paket: ‘nlme’
Das folgende Objekt ist maskiert ‘package:dplyr’:
collapse
Lade nötiges Paket: rpart
spatstat.core 2.3-0
Lade nötiges Paket: spatstat.linnet
spatstat.linnet 2.3-0
spatstat 2.2-0 (nickname: ‘That's not important right now’)
For an introduction to spatstat, type ‘beginner’
library(pracma)
Attache Paket: ‘pracma’
Das folgende Objekt ist maskiert ‘package:spatstat.core’:
rat
Die folgenden Objekte sind maskiert von ‘package:spatstat.geom’:
integral, quad
Das folgende Objekt ist maskiert ‘package:purrr’:
cross
library(plyr)
----------------------------------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------------------------------
Attache Paket: ‘plyr’
Die folgenden Objekte sind maskiert von ‘package:dplyr’:
arrange, count, desc, failwith, id, mutate, rename, summarise, summarize
Das folgende Objekt ist maskiert ‘package:purrr’:
compact
library(plotly)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Attache Paket: ‘plotly’
Die folgenden Objekte sind maskiert von ‘package:plyr’:
arrange, mutate, rename, summarise
Das folgende Objekt ist maskiert ‘package:ggplot2’:
last_plot
Das folgende Objekt ist maskiert ‘package:stats’:
filter
Das folgende Objekt ist maskiert ‘package:graphics’:
layout
library(mltools)
Attache Paket: ‘mltools’
Das folgende Objekt ist maskiert ‘package:tidyr’:
replace_na
library(data.table)
data.table 1.14.0 using 2 threads (see ?getDTthreads). Latest news: r-datatable.com
Attache Paket: ‘data.table’
Das folgende Objekt ist maskiert ‘package:spatstat.geom’:
shift
Die folgenden Objekte sind maskiert von ‘package:dplyr’:
between, first, last
Das folgende Objekt ist maskiert ‘package:purrr’:
transpose
library(htmlwidgets)
Function to calculate B1 and B2 of a specific omega.
calcB2 <- function(X1w, X2w, X3w) {
norm12 = norm(X1w-X2w,"2")
norm13 = norm(X1w-X3w, "2")
norm23 = norm(X2w-X3w, "2")
N = c(norm12, norm13, norm23)
ind <- which.max(N)
if(ind == 1){
B2w = X3w
}else if(ind == 2){
B2w = X2w
}
else{
B2w = X1w
}
return(list("B2w" = B2w, "ind" = ind))
}
Calculate ordered random variables for a specific omega.
calcOrder <- function(X1w, X2w, X3w, B2w, ind){
norm1 = norm(X1w-B2w, "2")
norm2 = norm(X2w-B2w, "2")
norm3 = norm(X3w-B2w, "2")
X3wP <- B2w
if (ind == 1){
if(norm2 >= norm3){
X1wP <- X2w
X2wP <- X3w
}else{
X1wP <- X3w
X2wP <- X2w
}
}else if(ind == 2){
if(norm1 >= norm3){
X1wP <- X1w
X2wP <- X3w
}else{
X1wP <- X3w
X2wP <- X1w
}
}else{
if(norm1 >= norm2){
X1wP <- X1w
X2wP <- X2w
}else{
X1wP <- X2w
X2wP <- X1w
}
}
return(list("X1" = X1wP, "X2" = X2wP, "X3" = X3wP))
}
Calculate F_i for a specific Omega
calcF <- function(R2, B2w, X1wP, X2wP, X3wP){
# calculate F1
# browser()
n <- length(R2$x)
F1ind = numeric()
b <- Norm(X1wP - B2w)
for(i in 1:n){
# browser()
x = as.numeric(R2[i,])
a <- Norm(x - B2w)
if(a >= b){
F1ind <- c(F1ind,i)
}
}
F1 <- R2[F1ind,]
#R2 <- R2[-F1ind,]
# calculate F2
F2ind = numeric()
b <- Norm(X2wP - B2w)
# print(b)
if(n!=0){
for(i in 1:n){
if(!(i %in% F1ind)){
x = as.numeric(R2[i,])
a <- Norm(x - B2w)
if(a >= b){
F2ind <- c(F2ind,i)
}
}
}
}
F2 <- R2[F2ind,]
#R2 <- R2[-F2ind,]
#calculate F3
#n <- length(R2$x)
F3ind = numeric()
b <- Norm(X3wP - B2w)
# print(b)
if(n != 0){
for(i in 1:n){
if(!(i %in% F1ind) && !(i %in% F2ind)){
x = as.numeric(R2[i,])
a <- Norm(x - B2w)
if(a >= b){
F3ind <- c(F3ind,i)
}
}
}
}
F3 <- R2[F3ind,]
#R2 <- R2[-F3ind,]
#calculate F4 (the rest) and F0 (the empty set)
F4 <- R2
return(list("F1ind" = F1ind, "F2ind" = F2ind, "F3ind" = F3ind))
}
Calculate Probability of one i which is part of C
calcProb <- function(sizeOmega, ind_x, ind, allF){
count = 0
for(i in 1:sizeOmega){
Fi_ind <- allF[ind,i]
# print(x)
# print(nrow(merge(x,Fi)))
Fi_ind = unlist(Fi_ind)
if(ind_x %in% Fi_ind){
count = count + 1
}
}
# print(count/sizeOmega)
return(count/sizeOmega)
}
Calculate C for one x (F_i (x) still missing)
funcC <- function(ind_x, sizeOmega, allF, empCDF){
res = 0
for(ind in 1:3){
# print(res)
res = res + calcProb(sizeOmega, ind_x, ind, allF)*(empCDF[ind_x]+(ind-1))
}
# print(res)
return(res/4)
}
calcAllC <- function(){
allF = replicate(sizeOmega, calcAllF())
n <- length(R2$x)
Cxy <- numeric(n)
dt <- data.table(x = rnorm(100,-5,5), y = rnorm(100,-5,5))
empCDF <- empirical_cdf(dt, ubounds=CJ(x = xR2, y = yR2))
empCDF <- empCDF$CDF
for(i in 1:n){
Cxy[i] <- funcC(i, sizeOmega, allF, empCDF)
}
final <- R2
final$Cxy <- Cxy
return(final)
}
Calculate all the X_(i), F_i for all omega
calcAllF <- function(){
# X1w = runif(2,lowerBound, upperBound)
# X2w = runif(2,lowerBound, upperBound)
# X3w = runif(2,lowerBound, upperBound)
X1w = rpois(2,4)
X2w = rpois(2,4)
X3w = rpois(2,4)
ret <- calcB2(X1w,X2w,X3w)
B2w <- ret$B2w
ind <- ret$ind
ret <- calcOrder(X1w, X2w, X3w, B2w, ind)
X1wP <- ret$X1
X2wP <- ret$X2
X3wP <- ret$X3
scriptF <- calcF(R2, B2w, X1wP, X2wP, X3wP)
return (scriptF)
}
Some testing